What Was Tested: Do crop yields increase from randomly placed GCOs globally? (T# indicated random GCO)
What This Allows Me To Do: Identify if yield-distance relationships (near analysis) are a function of evolutionary escape or other variables.
What Are The Model Assumptions: GLS models are used as they have less assumptions about error structure and heteroscadicity.
#Barley -T01
#Select Data and Rescale Distance
null<-read.csv(file="./Near_Null/GCO_Near_T01.csv", header=T)
null$rescale_ND <- null$NEAR_DIST/(1000*1000)
barley.null<-null %>% select(mean_barle,rescale_ND)
barley.null<-barley.null %>% filter(mean_barle > 0, na.rm=TRUE)
barley.null$logBarley<-log10(barley.null$mean_barle)
#Plot
ggplot(barley.null, aes(x=rescale_ND, y=logBarley)) +geom_point() +geom_smooth(method="lm") +theme_justin +xlab("Distance From GCO (1000's km)") +ylab ("Log10 Yield")
## `geom_smooth()` using formula 'y ~ x'
Barley Near Null Model, T01
### Model - Barley T01
# regular OLS no variance structure
barley.ols <- gls(logBarley ~ rescale_ND, data = barley.null)
# varFixed (variance changes linearly with X)
barley.fixed <- update(barley.ols, .~., weights = varFixed(~rescale_ND))
# varPower (variance changes as a power function with X)
barley.power <- update(barley.ols, . ~ ., weights = varPower(form = ~rescale_ND))
# varExp (variance changes as an exponential function of x)
barley.exp <- update(barley.ols, . ~., weights = varExp(form = ~ rescale_ND)) # error
# varConstPower (constant plus a power function of X (useful if X includes 0))
barley.ConstPower <- update(barley.ols, . ~., weights = varConstPower(form = ~ rescale_ND))
# compare all models by AIC
AIC(barley.ols, barley.fixed, barley.power, barley.ConstPower, barley.exp)
## df AIC
## barley.ols 3 19860.72
## barley.fixed 3 20649.60
## barley.power 4 19818.06
## barley.ConstPower 5 19818.25
## barley.exp 4 19818.39
#varPower
summary(barley.power)
## Generalized least squares fit by REML
## Model: logBarley ~ rescale_ND
## Data: barley.null
## AIC BIC logLik
## 19818.06 19847.35 -9905.029
##
## Variance function:
## Structure: Power of variance covariate
## Formula: ~rescale_ND
## Parameter estimates:
## power
## 0.101525
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 0.09030799 0.016095330 5.610820 0
## rescale_ND -0.01134541 0.001242831 -9.128682 0
##
## Correlation:
## (Intr)
## rescale_ND -0.939
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -4.2685673 -0.4483371 0.2915946 0.7273243 1.5829412
##
## Residual standard error: 0.4567186
## Degrees of freedom: 11197 total; 11195 residual
#Cassava -T01
#Select Data and Rescale Distance
null<-read.csv(file="./Near_Null/GCO_Near_T01.csv", header=T)
null$rescale_ND <- null$NEAR_DIST/(1000*1000)
cassava.null<-null %>% select(mean_cassa,rescale_ND)
cassava.null<-cassava.null %>% filter(mean_cassa > 0, na.rm=TRUE)
cassava.null$logCassava<-log10(cassava.null$mean_cassa)
#Plot
ggplot(cassava.null, aes(x=rescale_ND, y=logCassava)) +geom_point() +geom_smooth(method="lm") +theme_justin +xlab("Distance From GCO (1000's km)") +ylab ("Log10 Yield")
## `geom_smooth()` using formula 'y ~ x'
Cassava Near Null Model, T01
### Model - Cassava T01
# regular OLS no variance structure
cassava.ols <- gls(logCassava ~ rescale_ND, data = cassava.null)
# varFixed (variance changes linearly with X)
cassava.fixed <- update(cassava.ols, .~., weights = varFixed(~rescale_ND))
# varPower (variance changes as a power function with X)
cassava.power <- update(cassava.ols, . ~ ., weights = varPower(form = ~rescale_ND))
# varExp (variance changes as an exponential function of x)
cassava.exp <- update(cassava.ols, . ~., weights = varExp(form = ~ rescale_ND)) # error
# varConstPower (constant plus a power function of X (useful if X includes 0))
cassava.ConstPower <- update(cassava.ols, . ~., weights = varConstPower(form = ~ rescale_ND))
# compare all models by AIC
AIC(cassava.ols, cassava.fixed, cassava.power, cassava.ConstPower, cassava.exp)
## df AIC
## cassava.ols 3 9128.050
## cassava.fixed 3 10246.056
## cassava.power 4 9129.169
## cassava.ConstPower 5 9131.167
## cassava.exp 4 9127.686
#varExp
summary(cassava.exp)
## Generalized least squares fit by REML
## Model: logCassava ~ rescale_ND
## Data: cassava.null
## AIC BIC logLik
## 9127.686 9154.688 -4559.843
##
## Variance function:
## Structure: Exponential of variance covariate
## Formula: ~rescale_ND
## Parameter estimates:
## expon
## 0.002698391
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 0.6418420 0.015944886 40.25379 0
## rescale_ND 0.0080824 0.001226775 6.58837 0
##
## Correlation:
## (Intr)
## rescale_ND -0.92
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -5.3384659 -0.3436627 0.2918642 0.7277743 1.5659812
##
## Residual standard error: 0.4814093
## Degrees of freedom: 6316 total; 6314 residual
#Groundnut -T01
#Select Data and Rescale Distance
null<-read.csv(file="./Near_Null/GCO_Near_T01.csv", header=T)
null$rescale_ND <- null$NEAR_DIST/(1000*1000)
groundnut.null<-null %>% select(mean_groun,rescale_ND)
groundnut.null<-groundnut.null %>% filter(mean_groun > 0, na.rm=TRUE)
groundnut.null$logGroundnut<-log10(groundnut.null$mean_groun)
#Plot
ggplot(groundnut.null, aes(x=rescale_ND, y=logGroundnut)) +geom_point() +geom_smooth(method="lm") +theme_justin +xlab("Distance From GCO (1000's km)") +ylab ("Log10 Yield")
## `geom_smooth()` using formula 'y ~ x'
Groundnut Near Null Model, T01
### Model - Groundnut T01
# regular OLS no variance structure
groundnut.ols <- gls(logGroundnut ~ rescale_ND, data = groundnut.null)
# varFixed (variance changes linearly with X)
groundnut.fixed <- update(groundnut.ols, .~., weights = varFixed(~rescale_ND))
# varPower (variance changes as a power function with X)
groundnut.power <- update(groundnut.ols, . ~ ., weights = varPower(form = ~rescale_ND))
# varExp (variance changes as an exponential function of x)
groundnut.exp <- update(groundnut.ols, . ~., weights = varExp(form = ~ rescale_ND)) # error
# varConstPower (constant plus a power function of X (useful if X includes 0))
groundnut.ConstPower <- update(groundnut.ols, . ~., weights = varConstPower(form = ~ rescale_ND))
# compare all models by AIC
AIC(groundnut.ols, groundnut.fixed, groundnut.power, groundnut.ConstPower, groundnut.exp)
## df AIC
## groundnut.ols 3 13518.83
## groundnut.fixed 3 14568.52
## groundnut.power 4 13516.17
## groundnut.ConstPower 5 13518.17
## groundnut.exp 4 13520.29
#varPower
summary(groundnut.power)
## Generalized least squares fit by REML
## Model: logGroundnut ~ rescale_ND
## Data: groundnut.null
## AIC BIC logLik
## 13516.17 13544.27 -6754.084
##
## Variance function:
## Structure: Power of variance covariate
## Formula: ~rescale_ND
## Parameter estimates:
## power
## 0.0340043
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) -0.2628995 0.016587741 -15.849027 0
## rescale_ND 0.0050640 0.001237187 4.093157 0
##
## Correlation:
## (Intr)
## rescale_ND -0.933
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -4.4907423 -0.4119282 0.2553681 0.6928837 1.6701106
##
## Residual standard error: 0.5010125
## Degrees of freedom: 8322 total; 8320 residual
#Maize -T01
#Select Data and Rescale Distance
null<-read.csv(file="./Near_Null/GCO_Near_T01.csv", header=T)
null$rescale_ND <- null$NEAR_DIST/(1000*1000)
maize.null<-null %>% select(mean_maize,rescale_ND)
maize.null<-maize.null %>% filter(mean_maize > 0, na.rm=TRUE)
maize.null$logMaize<-log10(maize.null$mean_maize)
#Plot
ggplot(maize.null, aes(x=rescale_ND, y=logMaize)) +geom_point() +geom_smooth(method="lm") +theme_justin +xlab("Distance From GCO (1000's km)") +ylab ("Log10 Yield")
## `geom_smooth()` using formula 'y ~ x'
Maize Near Null Model, T01
### Model - Maize T01
# regular OLS no variance structure
maize.ols <- gls(logMaize ~ rescale_ND, data = maize.null)
# varFixed (variance changes linearly with X)
maize.fixed <- update(maize.ols, .~., weights = varFixed(~rescale_ND))
# varPower (variance changes as a power function with X)
maize.power <- update(maize.ols, . ~ ., weights = varPower(form = ~rescale_ND))
# varExp (variance changes as an exponential function of x) - Does not Converge
#maize.exp <- update(maize.ols, . ~., weights = varExp(form = ~ rescale_ND)) # error
# varConstPower (constant plus a power function of X (useful if X includes 0))
maize.ConstPower <- update(maize.ols, . ~., weights = varConstPower(form = ~ rescale_ND))
# compare all models by AIC
AIC(maize.ols, maize.fixed, maize.power, maize.ConstPower) #,maize.exp)
## df AIC
## maize.ols 3 24154.27
## maize.fixed 3 26364.94
## maize.power 4 24151.90
## maize.ConstPower 5 24153.91
#varPower
summary(maize.power)
## Generalized least squares fit by REML
## Model: logMaize ~ rescale_ND
## Data: maize.null
## AIC BIC logLik
## 24151.9 24182.03 -12071.95
##
## Variance function:
## Structure: Power of variance covariate
## Formula: ~rescale_ND
## Parameter estimates:
## power
## -0.02496171
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 0.3224096 0.013719180 23.50064 0
## rescale_ND -0.0156234 0.001062368 -14.70622 0
##
## Correlation:
## (Intr)
## rescale_ND -0.933
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -4.6035233 -0.4830642 0.2177102 0.6990680 2.0207779
##
## Residual standard error: 0.6155705
## Degrees of freedom: 13793 total; 13791 residual
#Millet -T01
#Select Data and Rescale Distance
null<-read.csv(file="./Near_Null/GCO_Near_T01.csv", header=T)
null$rescale_ND <- null$NEAR_DIST/(1000*1000)
millet.null<-null %>% select(mean_mille,rescale_ND)
millet.null<-millet.null %>% filter(mean_mille > 0, na.rm=TRUE)
millet.null$logMillet<-log10(millet.null$mean_mille)
#Plot
ggplot(millet.null, aes(x=rescale_ND, y=logMillet)) +geom_point() +geom_smooth(method="lm") +theme_justin +xlab("Distance From GCO (1000's km)") +ylab ("Log10 Yield")
## `geom_smooth()` using formula 'y ~ x'
Millet Near Null Model, T01
### Model - Millet T01
# regular OLS no variance structure
millet.ols <- gls(logMillet ~ rescale_ND, data = millet.null)
# varFixed (variance changes linearly with X)
millet.fixed <- update(millet.ols, .~., weights = varFixed(~rescale_ND))
# varPower (variance changes as a power function with X)
millet.power <- update(millet.ols, . ~ ., weights = varPower(form = ~rescale_ND))
# varExp (variance changes as an exponential function of x)
millet.exp <- update(millet.ols, . ~., weights = varExp(form = ~ rescale_ND)) # error
# varConstPower (constant plus a power function of X (useful if X includes 0))
millet.ConstPower <- update(millet.ols, . ~., weights = varConstPower(form = ~ rescale_ND))
# compare all models by AIC
AIC(millet.ols, millet.fixed, millet.power, millet.ConstPower, millet.exp)
## df AIC
## millet.ols 3 12808.62
## millet.fixed 3 13255.59
## millet.power 4 12805.30
## millet.ConstPower 5 12807.30
## millet.exp 4 12807.72
#varPower
summary(millet.power)
## Generalized least squares fit by REML
## Model: logMillet ~ rescale_ND
## Data: millet.null
## AIC BIC logLik
## 12805.3 12832.88 -6398.648
##
## Variance function:
## Structure: Power of variance covariate
## Formula: ~rescale_ND
## Parameter estimates:
## power
## 0.0553718
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) -0.16989423 0.026443817 -6.424724 0
## rescale_ND -0.01371859 0.001853118 -7.402976 0
##
## Correlation:
## (Intr)
## rescale_ND -0.966
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -4.8936737 -0.4107742 0.2698067 0.6705882 1.7464964
##
## Residual standard error: 0.5033066
## Degrees of freedom: 7299 total; 7297 residual
#Rapeseed -T01
#Select Data and Rescale Distance
null<-read.csv(file="./Near_Null/GCO_Near_T01.csv", header=T)
null$rescale_ND <- null$NEAR_DIST/(1000*1000)
rapeseed.null<-null %>% select(mean_rapes,rescale_ND)
rapeseed.null<-rapeseed.null %>% filter(mean_rapes > 0, na.rm=TRUE)
rapeseed.null$logRapeseed<-log10(rapeseed.null$mean_rapes)
#Plot
ggplot(rapeseed.null, aes(x=rescale_ND, y=logRapeseed)) +geom_point() +geom_smooth(method="lm") +theme_justin +xlab("Distance From GCO (1000's km)") +ylab ("Log10 Yield")
## `geom_smooth()` using formula 'y ~ x'
Rapeseed Near Null Model, T01
### Model - Rapeseed T01
# regular OLS no variance structure
rapeseed.ols <- gls(logRapeseed ~ rescale_ND, data = rapeseed.null)
# varFixed (variance changes linearly with X)
rapeseed.fixed <- update(rapeseed.ols, .~., weights = varFixed(~rescale_ND))
# varPower (variance changes as a power function with X)
rapeseed.power <- update(rapeseed.ols, . ~ ., weights = varPower(form = ~rescale_ND))
# varExp (variance changes as an exponential function of x)
rapeseed.exp <- update(rapeseed.ols, . ~., weights = varExp(form = ~ rescale_ND)) # error
# varConstPower (constant plus a power function of X (useful if X includes 0))
rapeseed.ConstPower <- update(rapeseed.ols, . ~., weights = varConstPower(form = ~ rescale_ND))
# compare all models by AIC
AIC(rapeseed.ols, rapeseed.fixed, rapeseed.power, rapeseed.ConstPower, rapeseed.exp)
## df AIC
## rapeseed.ols 3 14872.85
## rapeseed.fixed 3 15691.04
## rapeseed.power 4 14852.02
## rapeseed.ConstPower 5 14853.45
## rapeseed.exp 4 14853.28
#varPower
summary(rapeseed.power)
## Generalized least squares fit by REML
## Model: logRapeseed ~ rescale_ND
## Data: rapeseed.null
## AIC BIC logLik
## 14852.02 14880.55 -7422.01
##
## Variance function:
## Structure: Power of variance covariate
## Formula: ~rescale_ND
## Parameter estimates:
## power
## 0.07621067
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) -0.09282676 0.015800536 -5.874912 0
## rescale_ND -0.01535533 0.001219035 -12.596297 0
##
## Correlation:
## (Intr)
## rescale_ND -0.935
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -4.2768892 -0.4279157 0.3348087 0.6379341 1.5684906
##
## Residual standard error: 0.447725
## Degrees of freedom: 9259 total; 9257 residual
#Rice -T01
#Select Data and Rescale Distance
null$ric<-read.csv(file="./Near_Null/GCO_Near_T01.csv", header=T)
null$rescale_ND <- null$NEAR_DIST/(1000*1000)
rice.null<-null %>% select(mean_rice_,rescale_ND)
rice.null<-rice.null %>% filter(mean_rice_ > 0, na.rm=TRUE)
rice.null$logRice<-log10(rice.null$mean_rice_)
#Plot
ggplot(rice.null, aes(x=rescale_ND, y=logRice)) +geom_point() +geom_smooth(method="lm") +theme_justin +xlab("Distance From GCO (1000's km)") +ylab ("Log10 Yield")
## `geom_smooth()` using formula 'y ~ x'
Rice Near Null Model, T01
### Model - Rice T01
# regular OLS no variance structure
rice.ols <- gls(logRice ~ rescale_ND, data = rice.null)
# varFixed (variance changes linearly with X)
rice.fixed <- update(rice.ols, .~., weights = varFixed(~rescale_ND))
# varPower (variance changes as a power function with X)
rice.power <- update(rice.ols, . ~ ., weights = varPower(form = ~rescale_ND))
# varExp (variance changes as an exponential function of x)
rice.exp <- update(rice.ols, . ~., weights = varExp(form = ~ rescale_ND)) # error
# varConstPower (constant plus a power function of X (useful if X includes 0))
rice.ConstPower <- update(rice.ols, . ~., weights = varConstPower(form = ~ rescale_ND))
# compare all models by AIC
AIC(rice.ols, rice.fixed, rice.power, rice.ConstPower, rice.exp)
## df AIC
## rice.ols 3 13881.81
## rice.fixed 3 15470.57
## rice.power 4 13881.02
## rice.ConstPower 5 13883.02
## rice.exp 4 13880.58
#varPower
summary(rice.power)
## Generalized least squares fit by REML
## Model: logRice ~ rescale_ND
## Data: rice.null
## AIC BIC logLik
## 13881.02 13909.18 -6936.512
##
## Variance function:
## Structure: Power of variance covariate
## Formula: ~rescale_ND
## Parameter estimates:
## power
## -0.02412247
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) -0.02183335 0.017118768 -1.275404 0.2022
## rescale_ND 0.01065429 0.001261042 8.448799 0.0000
##
## Correlation:
## (Intr)
## rescale_ND -0.937
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -4.8420834 -0.4142408 0.2172051 0.6685366 1.6159786
##
## Residual standard error: 0.5835948
## Degrees of freedom: 8428 total; 8426 residual
#Rye -T01
#Select Data and Rescale Distance
null<-read.csv(file="./Near_Null/GCO_Near_T01.csv", header=T)
null$rescale_ND <- null$NEAR_DIST/(1000*1000)
rye.null<-null %>% select(mean_rye_Y,rescale_ND)
rye.null<-rye.null %>% filter(mean_rye_Y > 0, na.rm=TRUE)
rye.null$logRye<-log10(rye.null$mean_rye_Y)
#Plot
ggplot(rye.null, aes(x=rescale_ND, y=logRye)) +geom_point() +geom_smooth(method="lm") +theme_justin +xlab("Distance From GCO (1000's km)") +ylab ("Log10 Yield")
## `geom_smooth()` using formula 'y ~ x'
Rye Near Null Model, T01
### Model - Rye T01
# regular OLS no variance structure
rye.ols <- gls(logRye ~ rescale_ND, data = rye.null)
# varFixed (variance changes linearly with X)
rye.fixed <- update(rye.ols, .~., weights = varFixed(~rescale_ND))
# varPower (variance changes as a power function with X)
rye.power <- update(rye.ols, . ~ ., weights = varPower(form = ~rescale_ND))
# varExp (variance changes as an exponential function of x)
rye.exp <- update(rye.ols, . ~., weights = varExp(form = ~ rescale_ND)) # error
# varConstPower (constant plus a power function of X (useful if X includes 0))
rye.ConstPower <- update(rye.ols, . ~., weights = varConstPower(form = ~ rescale_ND))
# compare all models by AIC
AIC(rye.ols, rye.fixed, rye.power, rye.ConstPower, rye.exp)
## df AIC
## rye.ols 3 14778.26
## rye.fixed 3 15379.32
## rye.power 4 14756.76
## rye.ConstPower 5 14758.76
## rye.exp 4 14772.60
#varPower
summary(rye.power)
## Generalized least squares fit by REML
## Model: logRye ~ rescale_ND
## Data: rye.null
## AIC BIC logLik
## 14756.76 14784.91 -7374.381
##
## Variance function:
## Structure: Power of variance covariate
## Formula: ~rescale_ND
## Parameter estimates:
## power
## 0.08782291
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) -0.3955762 0.018380040 -21.52205 0
## rescale_ND 0.0175924 0.001446465 12.16231 0
##
## Correlation:
## (Intr)
## rescale_ND -0.939
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -4.7054619 -0.5228168 0.3049936 0.6535288 1.8361489
##
## Residual standard error: 0.4701497
## Degrees of freedom: 8400 total; 8398 residual
#Sorghum -T01
#Select Data and Rescale Distance
null<-read.csv(file="./Near_Null/GCO_Near_T01.csv", header=T)
null$rescale_ND <- null$NEAR_DIST/(1000*1000)
sorghum.null<-null %>% select(mean_sorgh,rescale_ND)
sorghum.null<-sorghum.null %>% filter(mean_sorgh > 0, na.rm=TRUE)
sorghum.null$logSorghum<-log10(sorghum.null$mean_sorgh)
#Plot
ggplot(sorghum.null, aes(x=rescale_ND, y=logSorghum)) +geom_point() +geom_smooth(method="lm") +theme_justin +xlab("Distance From GCO (1000's km)") +ylab ("Log10 Yield")
## `geom_smooth()` using formula 'y ~ x'
Sorghum Near Null Model, T01
### Model - Sorghum T01
# regular OLS no variance structure
sorghum.ols <- gls(logSorghum ~ rescale_ND, data = sorghum.null)
# varFixed (variance changes linearly with X)
sorghum.fixed <- update(sorghum.ols, .~., weights = varFixed(~rescale_ND))
# varPower (variance changes as a power function with X)
sorghum.power <- update(sorghum.ols, . ~ ., weights = varPower(form = ~rescale_ND))
# varExp (variance changes as an exponential function of x)
sorghum.exp <- update(sorghum.ols, . ~., weights = varExp(form = ~ rescale_ND)) # error
# varConstPower (constant plus a power function of X (useful if X includes 0))
sorghum.ConstPower <- update(sorghum.ols, . ~., weights = varConstPower(form = ~ rescale_ND))
# compare all models by AIC
AIC(sorghum.ols, sorghum.fixed, sorghum.power, sorghum.ConstPower, sorghum.exp)
## df AIC
## sorghum.ols 3 18645.01
## sorghum.fixed 3 18952.22
## sorghum.power 4 18477.67
## sorghum.ConstPower 5 18479.67
## sorghum.exp 4 18506.01
#varPower
summary(sorghum.power)
## Generalized least squares fit by REML
## Model: logSorghum ~ rescale_ND
## Data: sorghum.null
## AIC BIC logLik
## 18477.67 18506.6 -9234.833
##
## Variance function:
## Structure: Power of variance covariate
## Formula: ~rescale_ND
## Parameter estimates:
## power
## 0.1996702
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 0.24963025 0.015136174 16.49229 0
## rescale_ND -0.03253063 0.001207039 -26.95078 0
##
## Correlation:
## (Intr)
## rescale_ND -0.923
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -4.9640782 -0.4787471 0.2296422 0.7081050 1.8790790
##
## Residual standard error: 0.3669401
## Degrees of freedom: 10230 total; 10228 residual
#Soybean -T01
#Select Data and Rescale Distance
null<-read.csv(file="./Near_Null/GCO_Near_T01.csv", header=T)
null$rescale_ND <- null$NEAR_DIST/(1000*1000)
soybean.null<-null %>% select(mean_soybe,rescale_ND)
soybean.null<-soybean.null %>% filter(mean_soybe > 0, na.rm=TRUE)
soybean.null$logSoybean<-log10(soybean.null$mean_soybe)
#Plot
ggplot(soybean.null, aes(x=rescale_ND, y=logSoybean)) +geom_point() +geom_smooth(method="lm") +theme_justin +xlab("Distance From GCO (1000's km)") +ylab ("Log10 Yield")
## `geom_smooth()` using formula 'y ~ x'
Soybean Near Null Model, T01
### Model - Soybean T01
# regular OLS no variance structure
soybean.ols <- gls(logSoybean ~ rescale_ND, data = soybean.null)
# varFixed (variance changes linearly with X)
soybean.fixed <- update(soybean.ols, .~., weights = varFixed(~rescale_ND))
# varPower (variance changes as a power function with X)
soybean.power <- update(soybean.ols, . ~ ., weights = varPower(form = ~rescale_ND))
# varExp (variance changes as an exponential function of x)
soybean.exp <- update(soybean.ols, . ~., weights = varExp(form = ~ rescale_ND)) # error
# varConstPower (constant plus a power function of X (useful if X includes 0))
soybean.ConstPower <- update(soybean.ols, . ~., weights = varConstPower(form = ~ rescale_ND))
# compare all models by AIC
AIC(soybean.ols, soybean.fixed, soybean.power, soybean.ConstPower, soybean.exp)
## df AIC
## soybean.ols 3 17381.04
## soybean.fixed 3 18239.75
## soybean.power 4 17340.76
## soybean.ConstPower 5 17342.76
## soybean.exp 4 17344.47
#varPower
summary(soybean.power)
## Generalized least squares fit by REML
## Model: logSoybean ~ rescale_ND
## Data: soybean.null
## AIC BIC logLik
## 17340.76 17369.82 -8666.382
##
## Variance function:
## Structure: Power of variance covariate
## Formula: ~rescale_ND
## Parameter estimates:
## power
## 0.09661652
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 0.05650602 0.015105583 3.740738 2e-04
## rescale_ND -0.02214085 0.001149337 -19.264017 0e+00
##
## Correlation:
## (Intr)
## rescale_ND -0.935
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -5.3616282 -0.4015994 0.3470334 0.6821040 1.6531372
##
## Residual standard error: 0.4338111
## Degrees of freedom: 10549 total; 10547 residual
#Sunflower -T01
#Select Data and Rescale Distance
null<-read.csv(file="./Near_Null/GCO_Near_T01.csv", header=T)
null$rescale_ND <- null$NEAR_DIST/(1000*1000)
sunflower.null<-null %>% select(mean_sunfl,rescale_ND)
sunflower.null<-sunflower.null %>% filter(mean_sunfl > 0, na.rm=TRUE)
sunflower.null$logSunflower<-log10(sunflower.null$mean_sunfl)
#Plot
ggplot(sunflower.null, aes(x=rescale_ND, y=logSunflower)) +geom_point() +geom_smooth(method="lm") +theme_justin +xlab("Distance From GCO (1000's km)") +ylab ("Log10 Yield")
## `geom_smooth()` using formula 'y ~ x'
Sunflower Near Null Model, T01
### Model - Sunflower T01
# regular OLS no variance structure
sunflower.ols <- gls(logSunflower ~ rescale_ND, data = sunflower.null)
# varFixed (variance changes linearly with X)
sunflower.fixed <- update(sunflower.ols, .~., weights = varFixed(~rescale_ND))
# varPower (variance changes as a power function with X)
sunflower.power <- update(sunflower.ols, . ~ ., weights = varPower(form = ~rescale_ND))
# varExp (variance changes as an exponential function of x)
sunflower.exp <- update(sunflower.ols, . ~., weights = varExp(form = ~ rescale_ND)) # error
# varConstPower (constant plus a power function of X (useful if X includes 0))
sunflower.ConstPower <- update(sunflower.ols, . ~., weights = varConstPower(form = ~ rescale_ND))
# compare all models by AIC
AIC(sunflower.ols, sunflower.fixed, sunflower.power, sunflower.ConstPower, sunflower.exp)
## df AIC
## sunflower.ols 3 15781.89
## sunflower.fixed 3 16014.89
## sunflower.power 4 15567.05
## sunflower.ConstPower 5 15545.23
## sunflower.exp 4 15546.02
#varPower
summary(sunflower.power)
## Generalized least squares fit by REML
## Model: logSunflower ~ rescale_ND
## Data: sunflower.null
## AIC BIC logLik
## 15567.05 15595.6 -7779.525
##
## Variance function:
## Structure: Power of variance covariate
## Formula: ~rescale_ND
## Parameter estimates:
## power
## 0.2137333
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) -0.19299200 0.013591393 -14.19957 0
## rescale_ND -0.01446698 0.001141075 -12.67838 0
##
## Correlation:
## (Intr)
## rescale_ND -0.907
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -4.5243282 -0.4099035 0.2536765 0.7261783 1.6669870
##
## Residual standard error: 0.3364162
## Degrees of freedom: 9303 total; 9301 residual
#Wheat -T01
#Select Data and Rescale Distance
null<-read.csv(file="./Near_Null/GCO_Near_T01.csv", header=T)
null$rescale_ND <- null$NEAR_DIST/(1000*1000)
wheat.null<-null %>% select(mean_wheat,rescale_ND)
wheat.null<-wheat.null %>% filter(mean_wheat > 0, na.rm=TRUE)
wheat.null$logWheat<-log10(wheat.null$mean_wheat)
#Plot
ggplot(wheat.null, aes(x=rescale_ND, y=logWheat)) +geom_point() +geom_smooth(method="lm") +theme_justin +xlab("Distance From GCO (1000's km)") +ylab ("Log10 Yield")
## `geom_smooth()` using formula 'y ~ x'
Wheat Near Null Model, T01
### Model - Wheat T01
# regular OLS no variance structure
wheat.ols <- gls(logWheat ~ rescale_ND, data = wheat.null)
# varFixed (variance changes linearly with X)
wheat.fixed <- update(wheat.ols, .~., weights = varFixed(~rescale_ND))
# varPower (variance changes as a power function with X)
wheat.power <- update(wheat.ols, . ~ ., weights = varPower(form = ~rescale_ND))
# varExp (variance changes as an exponential function of x)
wheat.exp <- update(wheat.ols, . ~., weights = varExp(form = ~ rescale_ND)) # error
# varConstPower (constant plus a power function of X (useful if X includes 0))
wheat.ConstPower <- update(wheat.ols, . ~., weights = varConstPower(form = ~ rescale_ND))
# compare all models by AIC
AIC(wheat.ols, wheat.fixed, wheat.power, wheat.ConstPower, wheat.exp)
## df AIC
## wheat.ols 3 22286.50
## wheat.fixed 3 23485.16
## wheat.power 4 22266.17
## wheat.ConstPower 5 22267.78
## wheat.exp 4 22267.96
#varPower
summary(wheat.power)
## Generalized least squares fit by REML
## Model: logWheat ~ rescale_ND
## Data: wheat.null
## AIC BIC logLik
## 22266.17 22296 -11129.09
##
## Variance function:
## Structure: Power of variance covariate
## Formula: ~rescale_ND
## Parameter estimates:
## power
## 0.06416588
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 0.05869809 0.014548018 4.034783 1e-04
## rescale_ND -0.00866151 0.001127852 -7.679655 0e+00
##
## Correlation:
## (Intr)
## rescale_ND -0.937
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -4.7828632 -0.4473925 0.2572152 0.6835164 1.7102002
##
## Residual standard error: 0.4933548
## Degrees of freedom: 12812 total; 12810 residual